Method for analysis of pressure response in underground formations

ABSTRACT

A reservoir pressure in an underground formation surrounding a well is analyzed based on a direct measurement of the pressure at the wall of the well using the permeability of mud cake on the wall of the well in the region in which the pressure measurement is made; determining the thickness of mud cake on the well of the well; determining the hydrostatic pressure in the well in the region in which the pressure measurement is made; calculating a pressure decay index from the mud cake permeability and thickness, the hydrostatic pressure and the measured pressure; and using the pressure decay index to analyze the measured pressure to derive the reservoir pressure.

TECHNICAL FIELD

This invention relates to methods for analysing the pressure response in an underground formation, such as might be measured from a borehole passing through the formation. In particular, the methods apply to such methods for use when the formation pressure is influenced by the supercharging effect.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefits of priority from Application Number PCT/GB2005/001820, entitled “METHOD FOR ANALYSIS OF PRESSURE RESPONSE IN UNDERGROUND FORMATIONS,” filed under the PCT on May 10, 2005, which is commonly assigned to assignee of the present invention and hereby incorporated by reference in its entirety.

BACKGROUND ART

Formation pressure measurements made from wells play an important role in the management of reservoirs of underground fluids such as oil and gas. Because of their dynamic nature formation pressure measurements provide essential information on well productivity and dynamic reservoir description both in exploration and exploitation scenarios. Static pressure data can be used to compute formation fluid density and contacts. This can be important to determine reserves. Pressure transient data on the other hand can be important for estimating permeability and heterogeneity and average reservoir pressure.

Traditionally, pressure transient testing has taken the form of Drill Stem Testing (DST) or conventional well testing in which a well is put under test for a relatively long duration. While these can be excellent ways to meet test objectives, environmental and cost considerations do not allow use these techniques at all times. Wireline and LWD tools have been developed to make probe-based formation pressure measurements to address this issue.

Wireline and while-drilling formation testers counter many of the restrictions imposed by conventional well tests. While the theory of pressure transient analysis is applicable to data obtained by such formation tests, they require formulation to account for additional effects. Specifically, formation testers can be used during measurement while drilling. However, interpretation of the pressure data acquired in this dynamic environment can be challenging. One of the difficulties arises due to supercharging which results from mud filtrate invasion and changes significantly over the duration of drilling. This results in an increase in sandface pressure which is over and above the reservoir pressure. Therefore, any calculation of initial pressure and permeability must take into account the supercharging effect.

While drilling, the well bore pressure is normally maintained at a pressure substantially greater than the formation pressure by the use of drilling fluids to control production of formation fluids into the well bore (the drilling fluids or ‘muds’ are pumped through the wellbore and are also used for cuttings transport, cleaning of the drill bit and chemical stabilisation of the well). When a producing zone is penetrated, the wellbore sandface (the region of the wellbore wall in the producing zone) is exposed to mud pressure and filtrate immediately invades the near wellbore region. A mud cake is formed when drilling fluid flows into the formation and solids are deposited at the surface of the wellbore. This process is referred to as static filtration. As the mud cake grows it eventually stabilizes to a maximum thickness. This is as a result of the shearing action of the mud circulation as well as the mechanical action of the rotating drill pipe. This process is known as dynamic filtration. During these processes a pressure gradient is established in the formation.

A schematic description of the pressure profile with supercharging effect is shown in FIG. 1. The pressure in the wellbore near the surface of the mud cake is at hydrostatic pressure (pm) but drops rapidly across the mud cake (pa) and then gradually reduces in the formation, approaching formation (farfield) pressure (pI) some distance away from the wellbore. This near wellbore elevation in pressure above the farfield is known as the supercharging effect. From the above it is clear that if a pressure transient measurement were taken soon after drilling, any interpretation technique would have to take into account the effect of supercharged pressure.

Various techniques have been proposed to address the supercharging effect. Examples can be found in U.S. Pat. No. 5,602,334, PROETT, Mark, et al. FORMATION TESTING IN THE DYNAMIC DRILLING ENVIRONMENT. SPWLA 45th Annual Logging Symposium. June 6-92004., PROETT, Mark, et al. Formation Testing In the Dynamic Drilling Environment. IADC/SPE Drilling Conference. 2-4 Mar. 2004., PROETT, Mark, et al. Supercharge Pressure Compensation Using a New Wireline Testing Method and Newly Developed Early Time Spherical Flow Model. SPE 36524. 6-9 Oct. 1996, p. 329-342., GOODE, Peter, et al. Analytical models for a multiple probe formation tester. SPE 20737 September 1990., GOODE, Peter, et al. Influence of an invaded zone on a multiprobe formation tester. SPE Formation Evaluation. March 1996, p. 31-40.

This invention aims to provide a method of interpreting formation measurements that can account for the effect of supercharging.

DISCLOSURE OF THE INVENTION

One aspect of this invention provides a method of analysing a reservoir pressure in an underground formation surrounding a well, comprising:

-   -   determining the permeability of mud cake on the wall of the well         in the region in which the pressure measurement is made;     -   determining the thickness of mud cake on the well of the well in         the region in which the pressure measurement is made;     -   determining the hydrostatic pressure in the well in the region         in which the pressure measurement is made;     -   measuring the formation pressure at the wall of the well;     -   calculating a pressure decay index from the mud cake         permeability and thickness, the hydrostatic pressure and the         measured pressure; and     -   using the pressure decay index to analyse the measured pressure         to derive the reservoir pressure.

The pressure decay index can be calculated using the following relationship

$\beta = \frac{k_{m}\left( {p_{m} - p_{a}} \right)}{l_{m}{k\left( {p_{a} - p_{I}} \right)}}$

wherein

k_(m)=mud cake permeability

l_(m)=mud cake thickness

p_(m)=hydrostatic pressure

p_(a)=measured pressure

p_(I)=reservoir pressure

k=permeability.

Preferably, the method further comprises deriving at least one of horizontal permeability, vertical permeability and productivity index of the well in the region of the measurement.

A method according to the invention can comprise estimating at least one parameter and using non-linear regression to modify this estimate until the calculated or derived parameters result in correspondence with measured parameters.

Typical inputs to the analysis include a calculated invasion rate derived from mud cake properties, transient pressure computations from reservoir fluid and rock properties, formation pressure tester probe configuration parameters, pressure sampling rate and duration, and pressure transient data obtained from the pressure measurement.

With these inputs, the method preferably comprises determining a goodness of fit of pressure transient data.

In one embodiment of a method according to the invention, mud and mud cake properties are used to calculate invasion rate, and this invasion rate is applied, together with reservoir fluid properties, tester and formation configuration data and test data to a model with regression to provide reservoir pressure, permeability and productivity parameters.

The methods according to the invention can be applied to measurements made with wireline or while drilling formation tester tools.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a schematic diagram of a formation with supercharging;

FIG. 2 shows a schematic structure of the mathematical formulae underlying the embodiment of the invention;

FIG. 3 is a flow diagram of an interpretation workflow incorporating a method according to the invention;

FIG. 4 shows a comparison of pressure response with and without supercharging;

FIG. 5 is the pressure difference plot;

FIG. 6 shows the sensitivity to change in permeability;

FIG. 7 shows the sensitivity to change in initial reservoir pressure; and

FIG. 8 shows the sensitivity to change in beta factor (pressure decay index).

MODE(S) FOR CARRYING OUT THE INVENTION

This invention applies to measurements of formation pressure made using wireline pressure measurement tools, such as the MDT of Schlumberger, or more recent formation pressure while drilling (FPWD) tools. These are not described in detail here at their performance and properties are well known. These tools generally operate by applying a test probe against the wall of the wellbore (sandface) through any mud cake that might be present, and making pressure measurements and, optionally, taking samples of the formation fluid through the probe. Such measurements typically obtain data in the form of pressure and flow development over a period of time.

Data from such measurements are obtained digitally and are typically analysed by means of dedicated software applications to provide an indication of the formation properties around the well. The method of the present invention is based on a series of mathematical formulae that are discussed in more detail below. Variations may be made to these formulae while still retaining the essential methodology of the invention.

A structure for the derived mathematical formulae underlying one embodiment of the invention is presented in FIG. 2. In this case, the most important component is the pressure calculator that combines, by superposition, the formation test pressure response, filtrate invasion pressure response and diffusion of initial supercharged pressure. The filtrate invasion rate calculator computes the invasion rate that is used by the pressure calculator. In this case, the parameters to be computed from the formation pressure test are horizontal permeability, vertical permeability and undisturbed reservoir pressure. These parameters are used to estimate the productivity index of the well (PI). In addition, the initial pressure decay factor β is also determined.

The equations, presented in more detail below, assume a single-phase approximation; that is, the mobility of the mud filtrate and the reservoir fluid are similar. This assumption is approximately true in many cases. The effect of any divergence from this approximation can be neglected because the radius of mud filtrate invasion is generally much smaller than the radius of investigation, even for a wireline formation test.

A flow diagram of an interpretation workflow incorporating a method according to the invention is presented in FIG. 3.

The pressure calculator is a forward model. It computes pressure response as a function of time based on input parameters, some of which it is desired to compute in the first place. The pressure calculator is therefore used in a non-linear regression loop starting with the first estimates of the parameters of interest. A first estimate of horizontal permeability is obtained from logs taken while drilling or subsequently, vertical permeability is defaulted to 10 percent of horizontal permeability, and the initial reservoir pressure is considered to be hydrostatic pressure. A first estimate of the decay factor is taken using the method described below. The non-linear regression module is a standard, gradient-based algorithm tuned for pressure transient interpretation. The final outcome is the matched formation test pressure and the tuned parameters.

Knowing the reservoir pressure the productivity index of the well can be computed using standard industry methods.

In summary, data entry consists of the following:

a) Mud cake properties to calculate invasion rate.

b) Reservoir fluid and rock properties for the transient pressure computation. This includes initial estimates of output parameters.

c) Formation tester probe configuration.

d) Sampling rate and duration.

e) Corresponding pressure transient data.

The program outputs the following parameters:

a) Reservoir horizontal and vertical permeability

b) Initial reservoir pressure

c) Well productivity index

d) Goodness of fit of the pressure transient data.

The workflow can handle multiple probes. Therefore, both pre-test and vertical interference test can be analyzed. Outside the regression loop the pressure calculator is used for test design. Supercharging effect is generally prominent in low permeability reservoirs. FIG. 4 shows a comparison of pressure response with and without supercharging. FIG. 5 is the pressure difference plot. It is clear that the pressure profile is not only displaced but also has a different shape. This means that the permeability estimated by using the standard formation tester model would be different from that obtained by the proposed model, thus reinforcing the need to use the correct model.

The sensitivity to change in permeability, initial reservoir pressure and beta factor (pressure decay index) is presented in FIGS. 6, 7 and 8. Note that the clear separation of stabilized pressure during buildup in FIG. 6. This shows the effect of filtrate loss during buildup. In FIG. 7 the curvature of the pressure curve increases with the increase in difference between the initial probe pressure and the initial reservoir pressure. Thus, the larger the supercharging effect, the higher will be the discrepancy when using a standard model. In addition to determining the initial distribution of pressure, the decay factor also determines how fast the supercharging effect diffuses in the reservoir. This is demonstrated in FIG. 8. These figures also show that the pressure response is sensitive to the above-mentioned parameters and therefore can be resolved using non-linear regression.

The pressure calculator hooked to a standard non-linear regression routine is used to test the interpretation workflow. The observed test data used in our case is generated synthetically with a-priori knowledge of the reservoir parameters. Three cases are investigated. The probe pressure at the start of the test is fixed at 4100 psi. In the first case, the horizontal and vertical permeabilities are perturbed from the known values but the initial pressure is fixed. In the second case all the three parameters are perturbed. A comparison of the two cases suggests that while an increase in the number of unknowns adversely affects the quality of match, it is still good enough for all practical purposes. In the third case the initial pressure and filtrate invasion terms are disabled; that is, the match is obtained with a standard model used in formation testing. The match obtained is extremely poor, which is a clear demonstration of the need for specialized models. The actual values and the match obtained for the three cases are illustrated in the table below:

TABLE 1 Match with Match with Perm + Match with Actual Perm Initial Pressure sample-only Value unknown unknown model Horizontal 1 1.01 1.15 5.0 Permeability (md) Vertical 0.1 0.99 0.95 2.0 Permeability (md) Initial 4000 4000 4002 4100 Pressure

A summary of the essential physical problem to be modelled and the mathematical formulation underlying the formulae used in the methods of the invention is given below. It will be appreciated that deviation may be made from these formulae while still staying within the scope of the invention.

The solution in this case is obtained by the application of successive integral transforms to the governing equations and to the associated initial and boundary conditions. A brief exposition of the problem is given below.

The medium is bounded by the cylinder r=a and extends to ∞ in the direction of r positive. There are two cases in the z direction. In the first, z is unbounded; that is, (−∞<z<∞), and in the second, the formation is of thickness h; that is, (0<z<h) and has a no-flow boundary condition at the upper and lower boundary. An initial condition is superimposed, given by p(r,θ,z,t ₀)=(p _(a) −p _(I))e ^(−β() r−a)+p _(I)  (1)

The supercharged pressure, p_(a), is greater than the reservoir pressure, p_(I), due to the invasion of mud filtrate into the reservoir. Equation 1 shows an exponential decline of pressure from sandface to reservoir. When time t=t₀=0 (beginning of the test), r=a (at the sandface), p (a,θ,z,0)=p_(a) and as r→∞, p=p_(I). This decline might be represented by any arbitrary function. However, considering the nature of pressure diffusion, the exponential representation is accurate enough. The decay factor β determines the curvature of the pressure profile in the reservoir and depends on fluid and rock properties. It is possible to approximate this factor through actual reservoir simulation. Alternatively, a more simplistic but straightforward approach is to determine it from actual transient tests by non-linear regression. Since p_(a) is measured, in theory, if the mud filtrate invasion process can be rigorously modelled, it should be possible to compute p_(I) without having to impose an initial condition of the type given by Equation 1. For example, if the invasion history is known and the near well bore reservoir description is fairly accurate a reservoir simulator can be used to compute the supercharged pressure. However, this process is laborious and, often, without sufficient reliable data to validate the model. Hence, the focus of the methods of this invention is to be able to take advantage of the recorded pressure transient data. The pressure transient data influenced by the supercharged formation is interpreted to obtain reservoir parameters as well as the initial reservoir pressure (p_(I)), without having to resort to complex workflows. Instead of modelling how the pressure builds up to p_(a) from p_(I), an initial condition is imposed that is simple but can be history matched to the measured pressure transient data.

This approach is illustrated with reference to two cases discussed below

Case 1.0: The medium is bounded by the cylinder r=a and extends to ∞ in the direction of r positive and

${\left( {{- \infty} < z < \infty} \right) \cdot \frac{\partial_{p}\left( {a,\theta,z,t} \right)}{\partial_{r}}} = {{- \left( \frac{\mu}{k} \right)_{a}}q_{M}}$ The initial pressure situation is: p(r,θ,z,t ₀)=(p _(a) −p _(I))e ^(−β(r−α)) +p _(I) ; p(a,θ,z,t ₀)=p _(a) A continuous source at [a, 0, z₀] is introduced and the resulting pressure disturbance left to diffuse through a semi-infinite homogeneous porous medium.

The solution in Laplace space is given by

$\begin{matrix} {{\overset{\_}{p}\left( {a,0,z,s} \right)} = {{\frac{2{q(s)}{\mathbb{e}}^{- {st}_{0}}}{\pi^{3}\phi\; c_{t}a^{2}\sqrt{\eta_{z}}}\sum\limits_{m = 0}^{\infty}} \ni_{m}{{\int_{0}^{\infty}{\frac{{\mathbb{e}}^{{- {{z - z_{0}}}}\sqrt{\frac{{\eta_{r}\xi^{2}} + s}{\eta_{z}}}}}{{\xi\left( {{\eta_{r}\xi^{2}} + s} \right)}\left\{ {{J_{M}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}{\xi++}}\frac{4\eta_{r}}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}{q_{M}(s)}{\int_{0}^{\infty}{\frac{1}{{\xi\left( {{\eta_{r}\xi^{2}} + s} \right)}\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}{\xi++}}\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{V_{0}\left( {\xi\; a} \right)}{\left( {{\eta_{r}\xi^{2}} + s} \right)\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\xi}}}}}}} + \frac{p_{I}}{s}}}} & (2) \end{matrix}$

and in real time

$\begin{matrix} {{\overset{\_}{p}\left( {a,0,z,t} \right)} = {{\frac{2{U\left( {t - t_{0}} \right)}}{\pi^{3}\phi\; c_{t}a^{2}\sqrt{{\pi\eta}_{z}}}\sum\limits_{m = 0}^{\infty}} \ni_{m}}} & (3) \\ {\mspace{70mu}{\int_{0}^{\infty}{\int_{0}^{t - t_{0}}{\frac{{q\left( {t - t_{0} - \tau} \right)}{\mathbb{e}}^{{{- \eta_{r}}\xi^{2}\tau} - \frac{{({z - z_{0}})}^{2}}{4\eta_{z}\tau}}}{\xi\sqrt{\tau}\left\{ {{J_{m}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\tau}{\mathbb{d}{\xi++}}}}}\quad} & \; \\ {\mspace{65mu}{\frac{4\eta_{r}}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}(s){\int_{0}^{\infty}{\int_{0}^{t}{\frac{{q_{M}\left( {t - \tau} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\xi\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\tau}{\mathbb{d}{\xi++}}}}}}} & \; \\ {\mspace{205mu}{{\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{{V_{0}\left( {\xi\; a} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\xi}}}} + p_{I}}} & \; \end{matrix}$

For constant q equation (3) reduces to

$\begin{matrix} {{{{{{\overset{\_}{p}\left( {a,0,z,t} \right)} = {{\frac{{U\left( {t - t_{0}} \right)}q}{\pi^{3}\phi\; c_{t}a^{2}\sqrt{{\pi\eta}_{r}\eta_{z}}}\sum\limits_{m = 0}^{\infty}} \ni_{m}{\int_{0}^{\infty}{\frac{{\mathbb{e}}^{{- {\xi{({z - z_{0}})}}}\sqrt{\frac{\eta_{r}}{\eta_{z}}}}}{\xi^{2}\left\{ {{J_{m}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}} \times \times \left\{ {2 - {{\mathbb{e}}^{2\;{\xi{({z - z_{0}})}}\sqrt{\frac{\eta_{r}}{\eta_{z}}}}{{erfc}\left( {{\xi\sqrt{\eta_{r}\left( {t - t_{0}} \right)}} + \frac{\left( {z - z_{0}} \right)}{2\sqrt{\eta_{z}\left( {t - t_{0}} \right)}}} \right)}} + {{erfc}\left( {{\xi\sqrt{\eta_{r}\left( {t - t_{0}} \right)}} - \frac{\left( {z - z_{0}} \right)}{2\sqrt{\eta_{z}\left( {t - t_{0}} \right)}}} \right)}} \right\}{\mathbb{d}{\xi++}}\frac{4\eta_{r}}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}{\int_{0}^{\infty}{\int_{0}^{t}{\frac{{q_{M}\left( {t - \tau} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\xi\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\tau}{\mathbb{d}{\xi++}}}}}}}}}\quad}\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{{V_{0}\left( {\xi\; a} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\xi}}}} + p_{I}};{z > z_{0}}} & (4) \end{matrix}$

If z<z₀ we interchange z and z₀, where

$\eta_{r} = {{\frac{k_{r}}{\phi\; c_{t}\mu}\mspace{14mu}{and}\mspace{14mu}\eta_{z}} = \frac{k_{z}}{\phi\; c_{t}\mu}}$ ${U\left( {t - t_{0}} \right)} = \left\{ {\begin{matrix} 0 & {t < t_{0}} \\ 1 & {t \geq t_{0}} \end{matrix} \ni_{m}{= \left\{ {{\begin{matrix} \frac{1}{2} & {m = 0} \\ 1 & {{m = 1},2,3,\ldots} \end{matrix}{G_{v}\left( {\xi\; r} \right)}} = {{\left\{ {{{Y_{\upsilon}^{\prime}\left( {\xi\; a} \right)}{J_{\upsilon}\left( {\xi\; r} \right)}} - {{J_{\upsilon}^{!}\left( {\xi\; a} \right)}{Y_{\upsilon}\left( {\xi\; r} \right)}}} \right\}\mspace{14mu}{and}{V_{0}\left( {\xi\; a} \right)}} = {{\int_{a}^{\infty}{{re}^{{- \beta}\; r}{G_{0}\left( {\xi\; r} \right)}\ {\mathbb{d}r}{J_{\upsilon}(x)}}} = {{\sum\limits_{n = 0}^{\infty}{\frac{\begin{matrix} \left( {- 1} \right)^{n} \\ \left( \frac{x}{2} \right)^{\upsilon + {2n}} \end{matrix}}{\begin{matrix} {{n!}\Gamma} \\ \left( {\upsilon + n + 1} \right) \end{matrix}}\begin{matrix} {{{Bessel}\mspace{14mu}{function}\mspace{14mu}{of}\mspace{14mu}{the}}\mspace{11mu}} \\ {{first}\mspace{14mu}{kind}\mspace{14mu}{of}\mspace{14mu}{order}\mspace{14mu}\upsilon} \end{matrix}{Y_{\upsilon}(x)}}} = {\frac{{{J_{\upsilon}(x)}{\cos\left( {\upsilon\;\pi} \right)}} - {J_{- \upsilon}(x)}}{\sin\left( {\upsilon\;\pi} \right)}\begin{matrix} {{Bessel}\mspace{14mu}{function}\mspace{14mu}{of}\mspace{14mu}{the}} \\ {{second}\mspace{14mu}{kind}\mspace{14mu}{of}\mspace{14mu}{order}\mspace{14mu}\upsilon} \end{matrix}}}}}} \right.}} \right.$

Case 2.0: The medium is bounded by the cylinder r=a and extends to ∞ in the direction of r positive and

${\left( {0 < z < h} \right) \cdot \frac{\partial_{p}\left( {a,\theta,z,t} \right)}{\partial_{r}}} = {{- \left( \frac{\mu}{k} \right)_{a}}q_{M}}$ The initial pressure situation is p(r,θ,z,t ₀)=(p _(a) −p _(I))e ^(−β() r−a)+p _(I) ; p(a,θ,z,t ₀)=p _(a) A continuous source at [a, 0, z₀] is introduced and the resulting pressure disturbance left to diffuse through a semi-infinite homogeneous porous medium.

The solution in Laplace space is given by

$\begin{matrix} {{\overset{\_}{p}\left( {a,0,z,s} \right)} = {{\frac{4{q(s)}{\mathbb{e}}^{- {st}_{0}}}{\pi^{3}{ha}^{2}\phi\; c_{t}}\sum\limits_{m = 0}^{\infty}} \ni_{m}{\sum\limits_{n = 0}^{\infty}{\ni_{n}{{\cos\left( \frac{n\;\pi\; z_{0}}{h} \right)}{\cos\left( \frac{n\;\pi\; z}{h} \right)} \times \times}}}}} & (5) \\ {\mspace{45mu}{\int_{0}^{\infty}{\frac{1}{{\xi\left( {{\eta_{r}\xi^{2}} + {\eta_{z}\left( \frac{n\;\pi}{h} \right)}^{2} + s} \right)}\left\{ {{J_{m}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}{\xi++}}}}\quad} & \; \\ {\mspace{31mu}{\frac{4}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}{q_{M}(s)}{\int_{0}^{\infty}{\frac{1}{{\xi\left( {{\eta_{r}\xi^{2}} + s} \right)}\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}{\xi++}}}}}\quad} & \; \\ {\mspace{104mu}{{\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{V_{0}\left( {\xi\; a} \right)}{\left( {{\eta_{r}\xi^{2}} + s} \right)\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}{\mathbb{d}\xi}}}} + \frac{p_{I}}{s}}} & \; \end{matrix}$

and in real time

$\begin{matrix} {{\overset{\_}{p}\left( {a,0,z,t} \right)} = {{\frac{U\left( {t - t_{0}} \right)}{\pi^{3}{ha}^{2}\phi\; c_{t}}\sum\limits_{m = 0}^{\infty}} \ni_{m}{\int_{0}^{\infty}{\frac{1}{\xi\left\{ {{J_{m}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}{\int_{0}^{t - t_{0}}{q\left( {t -} \right.}}}}}} & (6) \\ {{{{{\left. {{{\mspace{124mu}\quad}t_{0}} - \tau} \right) \times \times}\quad}\quad}\quad}\left\lbrack {{\Theta_{3}\left\{ {\frac{\pi\left( {z - z_{0}} \right)}{2h},{\mathbb{e}}^{{- {(\frac{\pi}{h})}^{2}}\eta_{z}\tau}} \right\}} +} \right.} & \; \\ {{{{\left. \mspace{95mu}{\Theta_{3}\left( {\frac{\pi\left( {z + z_{0}} \right)}{2h},{\mathbb{e}}^{{- {(\frac{\pi}{h})}^{2}}\eta_{z}\tau}} \right\}} \right\rbrack{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}{\mathbb{d}\tau}{\mathbb{d}{\xi++}}}\quad}\quad}\quad} & \; \\ {{{\mspace{76mu}{\frac{4}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}{\int_{0}^{\infty}{\int_{0}^{t}{\frac{{q_{M}\left( {t - \tau} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\xi\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\tau}{\mathbb{d}{\xi++}}}}}}\quad}\quad}\quad} & \; \\ {\mspace{225mu}{{\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{{V_{0}\left( {\xi\; a} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}{\mathbb{d}\xi}}}} + p_{I}}} & \; \end{matrix}$

For constant q equation (6) reduces to

$\begin{matrix} {{\overset{\_}{p}\left( {a,0,z,t} \right)} = {{\frac{{U\left( {t - t_{0}} \right)}q}{\pi^{3}{ha}^{2}\phi\; c_{t}}\sum\limits_{m = 0}^{\infty}} \ni_{m}{{\int_{0}^{\infty}{\frac{1}{\xi\left\{ {{J_{m}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{m}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}{\int_{0}^{t - t_{0}}{\left\lbrack {{\Theta_{3}\left\{ {\frac{\pi\left( {z - z_{0}} \right)}{2h},{\mathbb{e}}^{{- {(\frac{\pi}{h})}^{2}}\eta_{z}\tau}} \right\}} + {\Theta_{3}\left( {\frac{\pi\left( {z + z_{0}} \right)}{2h},{\mathbb{e}}^{{- {(\frac{\pi}{h})}^{2}}\eta_{z}\tau}} \right\}}} \right\rbrack{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}{\mathbb{d}\tau}{\mathbb{d}{\xi++}}\frac{4}{\pi^{2}a}\left( \frac{\mu}{k} \right)_{M}{\int_{0}^{\infty}{\int_{0}^{t}{\frac{{q_{M}\left( {t - \tau} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\xi\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\tau}{\mathbb{d}{\xi++}}\frac{2\left( {p_{a} - p_{I}} \right){\mathbb{e}}^{\beta\alpha}}{\pi\; a}{\int_{0}^{\infty}{\frac{{V_{0}\left( {\xi\; a} \right)}{\mathbb{e}}^{{- \eta_{r}}\xi^{2}\tau}}{\left\{ {{J_{0}^{\prime\; 2}\left( {\xi\; a} \right)} + {Y_{0}^{\prime\; 2}\left( {\xi\; a} \right)}} \right\}}\ {\mathbb{d}\xi}}}}}}}}}} + p_{I}}}} & (7) \end{matrix}$

Where

${\Theta_{3}\left\{ {x,Q} \right\}} = {1 + {2{\sum\limits_{n = 1}^{\infty}{Q^{n^{2}}{{\cos\left( {2\;{nx}} \right)}\;\left\lbrack {{Q} < 1} \right\rbrack}}}}}$ Elliptic theta function of the third kind

The mud filtrate invasion, q_(M), could be modelled using any analytical function of time. However, for simplicity it is assumed that the mud cake is relatively thin compared to the well diameter and linear Darcy's law can be applied. This gives,

$\begin{matrix} {q_{M} = {\frac{2\pi\;{ahk}_{m}}{\mu\; l_{m}}\left( {p_{m} - p_{a}} \right)}} & (8) \end{matrix}$ where k_(m) and l_(m) are the permeability and thickness of the mud cake respectively and μ is the mud filtrate viscosity.

The pressure decay factor β describes the decay of the reservoir from the supercharged sandface pressure to the initial reservoir pressure. An initial estimate of β can be derived by imposing continuity of flow across sandface, which is

$\begin{matrix} {\beta = \frac{k_{m}\left( {p_{m} - p_{a}} \right)}{l_{m}{k\left( {p_{a} - p_{I}} \right)}}} & (9) \end{matrix}$

The effect of tool storage can be incorporated in the pressure solution by applying Duhamel's principle.

Nomenclature:

a wellbore radius, m.

c_(t) compressibility, P_(a) ⁻¹

φ porosity, fraction.

h layer thickness, m.

J_(v) vth order Bessel function of first kind.

J′_(v) derivative of vth order Bessel function of first kind.

Y_(v) vth order Bessel function of second kind.

Y′_(v) derivative of vth order Bessel function of second kind.

k_(m) mud filtrate permeability, m²

k_(r) horizontal permeability, m²

k_(z) vertical permeability, m²

l_(m) mud filtrate thickness, m.

μ viscosity, P_(a)·s.

p pressure, P_(a).

p_(a) probe pressure, P_(a).

p_(I) reservoir pressure, P_(a).

p_(m) hydrostatic pressure, P_(a).

q sampling rate, m³/s.

θ₃ elliptic theta function of the third kind.

q_(M) invasion rate, m³/s.

s Laplace variable.

t time, s.

t₀ time at start of test, s. 

The invention claimed is:
 1. A method for analysing a reservoir pressure in an underground formation surrounding a well, comprising: measuring a pressure at a wall of the well; determining a permeability of mud cake on the wall of the well in the region in which the pressure measurement is made; determining a thickness of mud cake on the well of the well in the region in which the pressure measurement is made; determining a hydrostatic pressure in the well in the region in which the pressure measurement is made; calculating a pressure decay index from the mud cake permeability and thickness, the hydrostatic pressure and the measured pressure; and using the pressure decay index and the measured pressure to derive the reservoir pressure.
 2. A method as claimed in claim 1, wherein the pressure decay index is calculated using the following relationship $\beta = \frac{k_{m}\left( {p_{m} - p_{a}} \right)}{l_{m}{k\left( {p_{a} - p_{I}} \right)}}$ wherein k_(m) =mud cake permeability l_(m) =mud cake thickness p_(m) =hydrostatic pressure p_(a) =measured pressure p_(I) =reservoir pressure k =permeability.
 3. A method as claimed in claim 1, further comprising deriving at least one of horizontal permeability, vertical permeability and a productivity index of the well in the region of the pressure measurement.
 4. A method as claimed in claim 1, further comprising using mud and mud cake properties to calculate an invasion rate, and applying this invasion rate together with reservoir fluid properties, tester and formation configuration data and test data to a model with regression to provide reservoir pressure, permeability and productivity parameters.
 5. A method as claimed in claim 1, wherein the pressure measurement comprises a sandface pressure measurement.
 6. A method as claimed in claim 1, wherein inputs to the analysis include a calculated invasion rate derived from mud cake properties, transient pressure computations from reservoir fluid and rock properties, formation pressure tester probe configuration parameters, pressure sampling rate and duration, and pressure transient data obtained from the pressure measurement.
 7. A method as claimed in claim 6, further comprising determining a goodness of fit of pressure transient data. 